clc
clear

% 原子半径
r0 = 0.55;
r1 = 0.45;

% 连接中心原子的原子或电子对个数个数. N_PAIRS=0,1,2,3,4
N_PAIRS=4;
% 其中孤对电子个数. N_ELECTRONS=0,1,...,N_PAIRS
N_ELECTRONS=2;

% 参考：w站词条“分子结构”
% 以下列出一些示例
% N_PAIRS, N_ELECTRONS, 分子形状, 代表性物质
% 2,0: 直线形，例如CO2
% 3,0: 三角形，例如BF3
% 4,0: 四面体形，例如CH4
% 4,1：三角锥形，例如NH3
% 4,2: 角形，例如H2O

function h = SPHERE(x,y,z,r)
  [X,Y,Z] = sphere(20);
  h = surf(r*X+x,r*Y+y,r*Z+z);
  set(h,'EdgeColor','none')
end

Atom = zeros(5,3);
isTransparent = zeros(5,1);

switch N_PAIRS
case 0
	Atom(1,:) = [0,0,0];
case 1
  a = r0+r1;
  Atom(1,:) = [0,0,0];
  Atom(2,:) = [a,0,0];
case 2
  a = r0+r1;
  Atom(1,:) = [0,0,0];
  Atom(2,:) = [a,0,0];
  Atom(3,:) = [-a,0,0];
case 3
  a = sqrt(3)*( r0+r1);
  Atom(1,:) = [a/2,sqrt(3)/6*a,0];
  Atom(2,:) = [0,0,0];
  Atom(3,:) = [a,0,0];
  Atom(4,:) = [a/2,sqrt(3)/2*a,0];
case 4
  a = sqrt(8/3)*(r0+r1);
  Atom(1,:) = [a/2,sqrt(3)/6*a,sqrt(6)/12*a];
  Atom(2,:) = [0,0,0];
  Atom(3,:) = [a,0,0];
  Atom(4,:) = [a/2,sqrt(3)/2*a,0];
  Atom(5,:) = [a/2,sqrt(3)/6*a,sqrt(6)/3*a];
otherwise
	error('invalid N_PAIRS');
	return;
end

if N_ELECTRONS > N_PAIRS
	error('invalid N_ELECTRONS.')
	return;
end
%
if N_ELECTRONS == 0
		;
else
	isTransparent(1+N_PAIRS-N_ELECTRONS+1:1+N_PAIRS) = 1;
end

% 平移到原点
for i = 2:N_PAIRS+1
  Atom(i,:) -= Atom(1,:);
end
Atom(1,:) = [0,0,0];

figure()
hold on
axis equal
for i = 1:N_PAIRS+1
  if i == 1
    r = r0;
  else
    r = r1;
  end
  h(i) = SPHERE(Atom(i,1),Atom(i,2),Atom(i,3),r);
  if isTransparent(i) == 1
    set(h(i),'FaceAlpha',0.3);
  end
end

if N_PAIRS >= 3
  for i = 2:N_PAIRS+1
    for j = i+1:N_PAIRS+1
      line([Atom(i,1) Atom(j,1)], [Atom(i,2) Atom(j,2)], [Atom(i,3) Atom(j,3)]);
    end
  end
end

